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Optimized 3D simulation method for modeling of 
out-of-plane radiation in silicon photonic integrated 

circuits 

Wouter J. Westerveld, H. Paul Urbach, and Mirvais Yousefi 


Abstract —We present an accurate and fast 3D simulation 
scheme for out-of-plane grating couplers, based on two dimen¬ 
sional rigorous (finite difference time domain) grating simu¬ 
lations, the effective index method (EIM), and the Rayleigh- 
Sommerfeld diffraction formula. In comparison with full 3D 
FDTD simulations, the rms difference in electric field is below 5% 
and the difference in power fiux is below 3%. A grating coupler 
for coupling from a silicon-on-insulator photonic integrated 
circuit to an optical fiber positioned 0.1 mm above the circuit is 
designed as example. 

Index Terms —Electromagnetic scattering by periodic struc¬ 
tures, FDTD methods, gratings, integrated optics, silicon on insu¬ 
lator technology, optical planar waveguide couplers, simulations. 


1. Introduction 

ILICON PHOTONICS became one of the most promising 
platform for integrated photonics, due to the possibil¬ 
ity of CMOS fabrication and thus the economy of scale, 
that came available in the last decade. While the original 
driver behind integrated optics was telecommunications, novel 
application fields such as on-chip interconnects and sensor 
systems have discovered the virtues of silicon photonics [1], 
[2]. Silicon photonic integrated circuits have a large refractive 
index contrast allowing for small device footprint, which is 
handy for most applications. However, it makes in-and-out 
coupling of light into the photonic integrated circuit (PIC) very 
difficult, since one has to match a ^9 jim fiber core with a 
~0.5 jim waveguide. The use of out-of-plane grating couplers, 
a technique already reported in the 1970’s, circumvent this 
problem by employing a grating that redirects the light from 
the waveguide upwards [3]. Radiation occurs from an area 
on the top surface of the PIC, allowing for a coupler with 
the same dimensions as the fiber. Alignment tolerances are 
large, (^7 jim) for an optical fiber positioned 102 jim above 
the grating. Grating couplers have seen a revival due to 
the ascendance of CMOS fabricated silicon photonics. The 
large index contrast allows for broadband coupling of light 
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(^25 nm). The gratings itself have also evolved over time from 
the original uniform ID gratings to today’s photonic crystal 
based 2D gratings [4], [5], [6], [7], [8]. 

The out-of-plane coupler is not only useful for coupling 
light into a fiber. Due to silicon’s lack of a direct band gap, 
generation or detection of light requires additional material. A 
near-future candidate as light source for mass-produced silicon 
PICs is a VCSEL mounted above an out-of-plane coupler with 
automated pick-and-place equipment [9], [10]. Furthermore, 
these couplers are very suitable for functional wafer-scale 
testing of PICs during the fabrication process. In the field 
of sensing, we expect grating couplers to be used for line-of- 
sight remote sensing in rough environments. This versatility of 
applications, in combination with the huge number of design 
parameters, requires a fast simulation method to calculate the 
full electromagnetic field at arbitrary position from the grating. 

Rigorous simulation methods such as finite difference time 
domain (FDTD) or eigenmode expansion (EMM) are required 
to accurately model the behavior of the high-contrast gratings. 
Current simulations of grating couplers are, to our knowledge, 
only performed as a 2D analysis (i.e. ID gratings). Two- 
dimensional gratings are understood as a superposition of 
two ID gratings [4], [5], [6], [7], [8]. In this paper, we 
show that this two-dimensional analysis is valid only in the 
vicinity of the grating coupler and we extend the method with 
three-dimensional free-space propagation. The focus of this 
work is mainly on uniform gratings to demonstrate our novel 
calculation technique, but the analysis is easily applicable to 
virtually any type of out-of-plane couplers, except for couplers 
with two-dimensional focusing of the light. 

The next section will detail the theory. In section III, we 
propose a novel simulation approach to describe the behavior 
of the grating coupler in three dimensions and compare it 
against full 3D FDTD simulations. In section IV, we use this 
method to design a uniform grating coupler within the limits 
of CMOS fabrication and to check the validity of full 2D 
simulations. 

H. Theory 

An out-of-plane grating coupler couples light from a high 
refractive index waveguide upwards via the air into, for 
example, an optical fiber. Figure 1 shows a schematic of an 
out-of-plane grating coupler, as well as the simulation scheme 
we propose in this paper. 
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(a) side-view schematics 


plane y2=7.3pm 



Fig. 1: (a) Schematics of the simulation approach, side-view 
of the grating coupler. The grating period A, the fill factor 
a and the etch depth b are indicated. The first three ’’point- 
sources” of the intuitive description are indicated P. The planes 
yi = 1.8 fim and y 2 = 7.3 ym are indicated. The numerical 
methods for the different steps of the simulation approach are 
given. 

(b) side-view and (c) top-view of the FDTD setup. E is the 
excitor. A1 and A2 are the absorbing shields. WG indicates the 
waveguide. GC indicates the grating coupler. SI is the sensor 
for incoming and back-reflected light through the waveguide. 
S2 is the sensor for transmission through the waveguide after 
the grating coupler. S3 is the sensor downward for radiation 
into the substrate. S4 is the sensor for upward radiation into 
the air. S5 is the sensor at the position of the fiber facet (full 
FDTD simulation). 


We consider the coupling of light from the fundamental 
mode of the input waveguide to the fundamental mode of 
the fiber. The simulation scheme we propose consists of four 
parts: a 2D FDTD simulation which describes the propagation 
from the waveguide to a plane just above the coupler. Then, 
the effective index method is applied to calculate the lateral 
profile of the field, based on the width of the grating in this 
plane, resulting in a 3D field profile. Thereafter, Rayleigh- 
Sommerfeld diffraction is used to propagate the field from 
this plane to the plane of the fiber facet and finally an overlap 
integral is used to calculate the coupling into the fiber mode. 
In what follows, we will detail the theory behind each of these 
steps. 

The behavior of the coupler can intuitively be understood 
by considering all the tall-short interfaces on the left edge of 
the grating grooves as ’’point-sources” which have a phase 


difference dictated by the propagation speed of the light 
through the waveguide (a few of such ’’point sources” are 
indicated by P in Fig. la). The effective (refractive) index of 
the grating can be estimated as the spatially weighted average 
of the effective indices of the fundamental modes in the tall 
and short parts of the waveguide. The fields emitted by these 
point-sources constructively interfere to a plane wave radiating 
under a certain angle Oq w.r.t. the y-axis. The relation between 
this angle and the waveguide properties is given by [11]: 

ns sin(6»g) = - q^, (1) 

where ns is the refractive index of the upper medium (air in 
Figure 1), neff is the effective index of the grating, q is the 
coupling order, Aq is the wavelength of the light in vacuum and 
A is the grating period. For perfect vertical coupling, Oq = 0°, 
Eq. (1) describes the second order resonance of a distributed 
Bragg reflector (DBR) [11], which very efficiently refiects the 
forward propagating light in the waveguide backwards, rather 
than radiating it upwards. Equation (1) can also be derived by 
treating the grating as a perturbation of the waveguide and 
calculating the coupling coefficient between the waveguide 
mode and a plane wave radiating with angle 6>q. The coupling 
strength scales with the fill factor a as sin(7ra) and depends 
linearly on the etch depth b [12]. For efficient coupling to 
a fiber, the grating should be designed such that only one 
coupling order (q = 1) is allowed, and the coupling strength 
should be designed between strong upwards coupling and 
low backwards refiections. The ratio between upwards and 
downwards radiation is strongly infiuenced by the buried oxide 
(BOX) height. For higher coupling efficiency, a non-uniform 
grating can be designed [5], [6]. 

A. Two-dimensional calculations and the effective index 
method 

The width of the vertical grating coupler is much larger 
than its height, thus the variation of the electromagnetic 
fields of the fundamental mode of the coupler is slow in 
the x-direction compared to its variation in the (y,z)-plane. 
A two-dimensional (y,z) analysis can therefore be used in the 
vicinity of the coupler up to approximately one wavelength 
above the grating. To obtain the lateral profile in a plane 
just above the coupler, we apply a method similar to the 
effective index method (EIM) by approximating the field E as 
E(x, y^ z) = E^^(^, z) • where E^^ is the electric 

field as calculated using two-dimensional analysis (assuming 
invariance in the x-direction) and lateral field profile 
is approximated from the lowest order mode in the waveguide. 


B. Propagation into the upper-half space 

The Rayleigh-Sommerfeld diffraction equation describes the 
electromagnetic field in a semi-infinite homogeneous medium, 
irradiated from a finite aperture [13]. In our case, the finite 
aperture is a sufficiently large part of the plane just above 
the grating coupler and air is the homogeneous medium. The 
field outside this aperture is assumed to be zero. The diffrac¬ 
tion equation for monochromatic light (vacuum wavelength 
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Ao = 27r//co) with free-space propagation constant k = k^n^ 
is written as [13]: 


U{x,y,z) = JJ U{x',yo,z')G{x,y,z;x',yo,z')dx'dz', 


aperture 


where 


( 2 ) 

(3) 

= J{x- x')^ + {y- y'oY + {z- z')"^, (4) 


277 r 


and U is any electric or magnetic field component in pha- 
sor notation with time dependence given by Greens’ 

function G is the sum of the fields of two in phase point 
sources that are images of each other with the plane y = y'o 
as mirror. This choice of Greens’ function allows the field in 
the air to be described as a function of the field in the aperture, 
not requiring knowledge of the divergence of the field. When 
the field is calculated for a horizontal plane (y constant), G 
depends on x — x' and 2 ; — z', so Eq. (2) is a 2D convolution, 
which can be very efficiently calculated using fast-fourier- 
transforms (FFTs) [14]. When the plane is tilted along the 
x-direction, the inner integral over x' is a convolution and 
can be calculated using FFTs. The 2D equivalent of Eq. (2), 
where U = U{y^z), is obtained by integrating Eq. (2) along 


x', resulting in [15]: 


U{y,z) = J U{yo,z')G{y,z-,yQ,z')dz', 

(5) 

aperture 


where 


G= 

(6) 

II 

1 

to 

+ 

1 

to 

(7) 


where H[^\kr) is the first-order Hankel function of the first 
kind, i.e. H[^\kr) = Ji{kr) -\-iYi{kr) where Ji and Yi are 
the first order Bessel functions of the first and second kind, 
respectively. 


C. Coupling into a fiber 

For monochromatic waves in dielectrics, the time-averaged 
power flux $ through a surface can be calculated as [16]: 

$ = 1/ Re(Ex-B*)dS, (8) 

^ J surface V k' J 

where E and B are the phasor notations of the electric and 
magnetic fields, respectively. The power coupling efficiency 
into a fiber mode with electric field (x, p) is estimated as 
[16]: 

_ \yEi{x,p)-Ei{x,p)ds\^ 

overlap \Ei.(x, p)\^dS ■ // \EI(x, p)\^dS' 

where E* is the electric field incident on the fiber facet. 
Coordinates x and p are in the tilted plane, parallel to the fiber 
facet. The fiber mode is approximated as a Gaussian beam 
with a beam diameter 2wq = 10.4 pm. Details of coupling 
from air into the fiber are neglected since there is a small 
refractive index step and a small angle of incidence of the 


TABFE I: Simulation speeds and memory requirements for 
the different numerical methods as presented in Sec. III. All 
simulations are performed with a 44 nm numerical gridsize. 
Note that the simulations in Sec. IV have a numerical gridsize 
of 20 nm or 10 nm, when necessary. Rayleigh-Sommerfeld 
(R.-S.) calculations are implemented in Matlab. This test 
is performed on 32bit PC with Windows XP professional, an 
Intel Xeon CPU 5130 @ 2.00 GHz and 3.00GB of RAM. 


Method 

Window [ym] 

Time 

Memory 

3D FDTD (upto y 2 ) 

32.56 X 12.10 X II 

17h 

2.62GB 

3D FDTD (upto yi) 

32.56 X 12.10 X 7.7 

12h 

2.15GB 

2D FDTD (upto ^ 2 ) 

32.56 X 11 

2m 

80MB 

2D FDTD (upto yi) 

32.56 X 7 

Im 

20MB 

3D R.-S. (horizontal plane) 

surface 21.42 x 10.34 

2s 

< 100MB 

3D R.-S. (tilted plane) 

surface 21.42 x 10.34 

10m 

< 100MB 

2D R.-S. 

length 21.42 

Is 

<1MB 

Effective Index Method 

surface 21.42 x 10.34 

Is 

<11MB 



Fig. 2: Flux versus wavelength. Comparison of 3D FDTD 
simulations and 2D FDTD simulations. 


incoming wave. When the field El. can be separated in x and p 
dependence, i.e El.{x^ p) = {x) • EI^p[ p), then the overlap 

^overlap Can also be separated, i.e. //overlap = //overlap,a: ' //overlap,p- 

III. Simulation scheme 

We suggest a novel combination of numerical methods, 
see Fig. 1, that employs the theory described in the previous 
section. The electromagnetic field up to plane just above the 
coupler is calculated using 2D FDTD simulations and the 
lateral profile is obtained using the effective index method 
(EIM). Propagation from this plane upwards is performed 
using the 3D Rayleigh-Sommerf eld diffraction formula. The 
method, as well as individual steps, are compared against full 
3D FDTD simulations. Table 1 lists the calculation speeds of 
the different steps detailed in this section. 

The grating coupler we investigate in this section has a 
silicon-on-insulator (SOI) layer stack which consists of a 
silicon substrate, covered by a 1980 nm thick Si02 layer for 
optical isolation between the substrate the 220 nm top mono¬ 
crystalline silicon device layer. The grating is partially etched 
into this device layer (see Fig. lb). The coupler waveguide 
has width W = 9.9 pm in the x-direction, the grating has 
period A = 616 nm, etch depth 6 = 88 nm, fill factor 
a = 50% and = 20 fundamental periods. The first groove 
of the grating is at z = 17.6 pm. FDTD simulations are 
performed in CrystalWave [17] with a 44 nm gridsize, which 
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Novel simulation scheme 





IT ■ II 



0 - 


Novel simulation scheme w.r.t. to 3D FDTD 

(d) Ex norm difference [V/|im] 




1 




^ -0.15 -0.1 -0.05 

0.05 

0.1 

0.15 


(e) Ex phase difference [rad] 



Fig. 3: Effective index method (EIM) compared to full 3D 
EDTD simulations. The plots present the Ex field in the bi¬ 
plane. Plots (a-c): Ex field as calculated using the EIM. Plots 
(d-f): difference between EIM and 3D EDTD simulations, 
i.e. E^^^ — Plots (a,b,d,e) show a top-view with 

the z-direction and the x-direction along the horizontal and 
vertical axes, respectively. The position of the grating coupler 
is indicated in the bottom of each plot. Plots (c and f) show 
a cross-section through x = 0. 


is a bit large for this grating but does describe all the relevant 
phenomena. A mode excitor launches a light pulse with central 
wavelength Aq = 1.55 yum and bandwidth 0.2 yum, forward 
traveling from z = 0.5 yum in the fundamental TE mode of 
the waveguide. In order to avoid artifacts of the excitor, we 
introduce two strongly absorbing and reflecting shields (A1 
and A2 in Eig. lb) in the system, such that at the point 
indicated WG in Eig. lb, all the power resides in the input 
waveguide. The power flux is normalized to the power flux 
incident on the grating coupler, as measured by sensor 1 (SI). 

The X component of the electric field. Ex, is much larger 
than the other components (3 orders of magnitude) and its 
real part is therefore used to measure the difference between 
alternative simulation methods (Ex) and a EDTD simulations 
(E^^^). Both phase and amplitude are taken into account. 


Eig. 4: Our novel simulation approach compared to full 3D 
EDTD simulations. The plots present the Ex field in the y 2 - 
plane. Eigure layout is identical to Eig. 3 


Eor the field in a surface, the root-mean-square difference is 
used as follows: 


^Ex 


ff Re{Ex - E^^^pdxdz 

surface 




ff Re{E^^^}‘^dxdz 

surface 


( 10 ) 


Another measure used is the relative difference in power flux 
4>. The latter takes all electric and magnetic field components 
into account. 

1) The effective index method: Two-dimensional simula¬ 
tions in the (y,z)-plane are used in the vicinity of the coupler 
up to a plane at = 1.8 yum above the grating, see Eig. 1 (b) 
and (c). The power flux through this plane is shown versus 
wavelength in Eig. 2 for both 2D and 3D simulations. The 
DBR behavior of the coupler is clearly visible at A = 1.61 yum 
as expected from Eq. (1). The difference between 2D and 3D 
simulations at central wavelength of 1.55 yum is below 1%. 

To obtain the lateral profile in the plane just above the 
coupler, we apply a method similar to the effective index 
method (EIM). The lateral field profile E^^^^(x) is approxi¬ 
mated from the lowest order mode in the waveguide which is 
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Fig. 5: Power coupling between the fundamental waveguide¬ 
mode and the fiber-mode. The distance of the PIC to the 
fiber (e.g. the fiber y-position) is plotted along the horizontal 
axis. For each y-position, the z-position and tilt of the fiber 
are optimized for maximum coupling. Grating parameters as 
specified in Sec. III. This result is calculated using our novel 
method presented in Sec. III. 


calculated using the film mode matching (FMM) mode solver 
as implemented in the FimmWave [17] software. Figure 3 
shows that the difference between the EIM (which uses 2D 
FDTD) and the full 3D FDTD, is = 5%. 

2) Rayleigh-Sommerfeld diffraction formula: Using the 
field at yi as obtained from a 3D FDTD simulation, we 
calculate the field at ^2 = 7.3 jim using the Rayleigh- 
Sommerfeld diffraction formula. We find that the difference in 
field as compared to a full 3D FDTD simulation, = 3%, 
and the difference in power flux, AT> = 2%. For full two- 
dimensional calculations using the 2D Rayleigh-Sommerfeld 
diffraction formula, Eq. (5), the differences with FDTD are 
AE^ = 2% and A^ = 1%. 

3) Novel simulation scheme: The results of our simulation 
approach are presented in Fig. 4. The difference in electric 
field is /S.Ex = 5% and the difference in the fiux is A4> = 3%. 
We therefore conclude that this approach is accurate enough to 
use in optimizing the design of grating couplers which couples 
light to fibers or to other optical components at arbitrary 
distance from the grating coupler. 

IV. Results 

This section presents results we obtained with our 3D 
simulations of out-of-plane grating couplers. In Sec. IV-A, we 
investigate up to which height the full 2D theory is valid. In 
Sec. IV-B, as a demonstration of our method, we design a 
grating which is optimized for coupling to an optical fiber 
that is positioned 102 ym above the chip. 

A. 2D versus 3D simulations 

Fully two-dimensional simulation approaches can be used 
to calculate fields in the vicinity of the grating coupler. The 
coupling efficiency is then calculated by applying the effective 
index method on the fiber facet, which allows to correct the 
2D simulation with a simple overlap efficiency in the lateral 
direction, r/overiap,^:^ as shown in Sec. II-C. Further away from 
the grating, the field will spread in lateral direction and a 
3D approach is required. In Fig. 5, we show the difference 
between 3D and 2D simulations. At each fiber height, the fiber 


(a) Power flux into air [dB] (b) Power flux into substrate [dB] 



500 550 600 650 700 500 550 600 650 700 

grating period [nm] grating period [nm] 


Fig. 6: Optimizing of the grating parameters, (a) Power 
radiated into the air. (b) Power radiated into the substrate, 
(c) Fiber tilt for maximum power coupling to the fiber, (d) 
Power coupling to the fiber. Result is calculated using our 
novel simulation approach. The grating parameters (A and h) 
for maximum upwards coupling and for maximum coupling to 
the fiber are indicated X. The dashed line indicates the grating 
parameters where a Bragg reflection is expected. The grating 
has a fill factor a = 50% and width W = 22 jim. 


z-position (longitudinal position) and tilt 0 are optimized for 
maximum coupling. The difference in to-fiber power coupling 
between a fully 2D analysis and a 3D analysis is ^2 dB when 
the fiber is 100 /im above the grating. 

B. Grating coupler design 

This section presents design steps for an out-of-plane 
grating coupler which couples light to a single mode fiber 
positioned 102 pm above the chip. In order to allow for 
mass-fabrication, we design a grating coupler that obey’s 
the latest design rules of industrial CMOS lines. The EU- 
FP7 consortium ePIXfab, see [18], offers access to high-end 
photonic IC technologies at IMEC and LETI and we will use 
their design rules [19], which require a 2 pm thick BOX layer 
and a minimal feature (groove or line) size of 120 nm. The 
height of our waveguide is 220 nm to suppress TM-modes. 
We found that the most dominant parameters in the grating 
design are the grating period A and the etch depth b. We will 
therefore start by optimizing the grating using these parameters 
and later fine-tune the system suing the other parameters as 
mentioned in Sec. I. 

The EDTD simulations have a numerical gridsize of 20 nm 
or, when required, 10 nm and is chosen such that it always 
perfectly represents the material profile of the grating coupler. 
Eree-space propagation has a numerical gridsize of 40 nm. The 
absorbing/refiecting shield is omitted as the numerical artifact 
of direct air radiation from the excitor does not influence 
the field at the fiber facet using this simulation scheme. Eor 
each set of grating parameters, the fiber position and tilt are 
optimized using a fully 2D approach. The first guess of the 
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Fig. 7: Power coupling between the fundamental waveguide 
mode and the fiber mode. The width of the waveguide is 
plotted along the horizontal axis. 



Fig. 8: Bandwidth of the grating coupler with parameters 
(A = 580 nm, 6 = 40 nm, a = 50%). 


fiber tilt for this optimization is given by Eq. (2) and the first 
guess of fiber z-position is the maximum Ex field component 
in the plane 102 /im above the coupler. The optimization steps 
are 100 nm in position and 0.1° in tilt. After finding the 
optimal fiber position, the coupling efficiency is calculated 
using a 3D simulation. 

A parameter scan of A and b in steps of 20 nm is presented 
in Fig. 6, with the fill factor a = 50% and the width W = 22 /j,m 
kept constant. The grating always has = 60 periods because 
there is only very little light left in the waveguide after this 
number of periods due to the strong coupling. Figures 6 (a) 
and (b) show the amount of upwards and downwards coupled 
light, which is calculated using the 2D FDTD simulation (also 
see Fig. 1). They show the same trend in dependence on 
A and 6, which makes it impossible to discriminate against 
downwards coupling without extra tricks. The dotted line is 
the position where the grating is expected to behave as a Bragg 
refiector, from Eq. (1) with 0i =0°. The dip in both up- and 
downwards fiux is clearly observed. Figure 6c shows the fiber 
tilt for optimal coupling into the fiber, this tilt is equal to the 
angle of the dominant plane-wave in the radiated field because 
the overlap of the fiber mode with the incident field is highest 
for a constant phase front. The prediction of the DBR line 
accurately describes vertical coupling. The difference between 
the optimal fiber tilt and the out-coupling angle 0i as predicted 
by Eq. (1) was below 3.5° for all simulations in this plot. The 
total coupling efficiency from the waveguide to the fiber is 
presented in Figure 6d. One sees that the grating parameters 
for maximum upwards coupling and maximum coupling to 
the fiber are not the same. In Fig. 6a, the maximal coupling 
is observed for period A = 580 nm and etch depth 6 = 40 nm. 
From this parameters, we optimized the duty cycle a, which 
was varied from 120 nm to 480 nm in steps of 20 nm. The 



-6 -4 -2 0 2 4 6 


fiber z offset [ |a.m] 

Fig. 9: Power coupling from the waveguide to the fiber for fiber 
alignment offsets as indicated by the axes. Grating coupler 
with parameters (A = 580 nm, 6 = 40 nm, a = 50%). 


coupling efficiency to the fiber showed only a weak change 
when varying a, except when the effective index of the grating 
was such that the DBR refiection showed up. The maximum 
was found at a = 50%. 

The next parameter being optimized is the width of the 
grating coupler waveguide, as presented in Fig. 7. The optimal 
width is found to be 22 fim. 

For this coupler, the bandwidth for a maximum transmission 
loss of 1 dB is 25 nm, as shown in Fig. 8. The fiber is 
102 fim above the coupler, so light of a different wavelength 
will radiate under a different angle, thereby not fully aiming 
at the fiber. The fiber position was given offsets in z- and x- 
directions to test the alignment tolerances. The area with a 
maximum transmission loss of IdB has a diameter of 6.8 //m 
and is almost circular (see Fig. 9). 

The optimal grating parameters we found are slightly differ¬ 
ent from the ones in Ref. [6]. This can be attributed to the fact 
that we chose the fiber to be positioned further away from the 
grating. Our etch depth is shallower, resulting in less strong 
coupling per propagated length through the coupler, and thus 
a larger radiative area. This gives rise to less divergence of the 
electromagnetic field when propagating upwards. To illustrate 
this. Figures 10 and 11 present the field at the plane of the 
fiber facet. The first figure shows the plane as irradiated by 
a grating coupler with a shallow etch depth, b = 40 nm. The 
second grating coupler as used for the second figure has a 
deeper etch depth, b = 10 nm. In both cases, the position of 
the fiber is optimized for maximum coupling, as in Fig. 6. 
The figures show that the field is slightly less diverged for 
the shallow grating and also that the phase of the field is 
much more constant over the positions where the amplitude 
of the field is strongest. The novel simulation method gives 
the opportunity to inspect these fields and do this analysis. 
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Fig. 10: Electric field (x component) at the plane of the fiber 
facet, (a) and (b) are a top-view with the p and the x along 
the horizontal and vertical axis, respectively. Plot (a) shows 
the amplitude of Ex. Plot (b) shows the phase of Ex. The 
overall phase of the field is rotated for clarity of the plot. 
Plot (c) is a cross-section of the field at x = 0 pm. The solid 
and dashed line are the amplitude (left axis) and phase (right 
axis), respectively. The field originates from a grating with 
parameters (A = 580 nm, 6 = 40 nm, a = 50%). 

V. Conclusion 

In this paper, we presented a novel three-dimensional sim¬ 
ulation scheme and used it to simulate out-of-plane grating 
couplers. The method is fast and accurate enough to meet 
the requirements for grating coupler design, involving a huge 
number of design parameters. In comparison with 3D FDTD 
simulations, our method showed a root-mean-square difference 
in electric field below 5% and a difference in power fiux below 
3%, while improving computational speed by two orders of 
magnitude. 

We foresee a broad range of applications for grating cou¬ 
plers in silicon photonics, and our method fulfills this versatile 
need. The presented method is applicable to may out-of-plane 
radiation simulations, such as coupling to the emission of a 
VCSEL, or on-wafer inspection of photonic integrated circuits. 

Current simulations of grating couplers are, to our knowl¬ 
edge, only performed as a 2D analysis. We found that this 
approximation is only valid in the vicinity of the coupler by 
comparing 2D and 3D FDTD simulations. 

As an example of our method, we designed a CMOS com¬ 
patible grating coupler for coupling of light from a waveguide 
to a fiber at a distance of 0.1 mm from the photonic integrated 
circuit. 

In this work, a FDTD method is used to calculate the 
electromagnetic fields in the vicinity of the grating coupler, 
however, this is not a limitation. Other methods such as the 
eigemode expansion technique (EME) can also be used in 
the presented scheme [4]. In the field of sensing, we expect 
line-of-sight remote readout of silicon photonic integrated 
circuits to become possible using grating couplers, opening the 
possibility of sensing in extreme environments without need 
for a physical interconnect such as a fiber. 


Eig. 11: Eigure layout is identical to Eig. 10. This field 
originates from a grating with parameters (A = 630 nm, 
h = 10 nm, a = 50%). 
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